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Least Action Principle for the Real-Time Density Matrix 
Renormalization Group 
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A kind of least action principle is introduced for the discrete time evolution of one-dimensional 
quantum lattice models. Based on this principle, we obtain an optimal condition for the matrix 
product states on succeeding time slices generated by the real-time density matrix renormaliza- 
tion group method. This optimization can also be applied to classical simulations of quantum 
circuits. We discuss the time reversal symmetry in the fully optimized MPS. 
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§1. Introduction 

The density matrix renormalization group (DMRG) 
method has been widely applied to calculations of eigen- 
states of low-dimensional quantum systems. 1-3 ' The 
method can be regarded as a numerical variational 
method, which optimizes position dependent matrix 
product state (MPS) by way of iterative improvements 
of local matrices. 4 ' 5 - ) This variational background guar- 
antees that truncation error in the block spin transfor- 
mation does not accumulate in the iterative numerical 
calculations of the finite system DMRG algorithm. 

One of the recent progress in DMRG is its applica- 
tion to quantum states under real or imaginary time 
evolution. 6 ~ 10 ) The concept of adopted time dependence 
is a key feature in the real-time DMRG method, where 
the evolving quantum state is approximated by MPS as 
precise as possible at each time slice. 8 ~ 11 - ) In this arti- 
cle we focus on the weak breaking of the time-reversal 
symmetry in the numerical algorithm of the real-time 
DMRG method. Suppose that we start from an ini- 
tial state \9(tj)). After the numerical time evolution 
with respect to the time-independent Hamiltonian H, 
we get to the calculated final state \^(t F )) that approx- 
imates exp[(t F — t I )H/ih]\if>(t I )). The backward nu- 
merical time evolution from \^{t F )) toward past direc- 
tion gives the calculated state \ty'(tj)) that approximates 
exp[(tj - t F )H/ih]\y(t F )). The state |*'(^)) thus ob- 
tained is not actually the same as the initial state \^(tj)). 
This is an example of the slight asymmetry in time in 
the real-time DMRG method. Qualitatively speaking, 
the discrepancy between ^(tj)) and |^'(^)) can be at- 
tributed to the accumulation of truncation error caused 
by the repeated use of time adopted renormalization pro- 
cesses. 

In order to recover the time-reversal symmetry, we in- 
troduce a kind of least action principle, which is related 
to MPSs on all the time slices. Also we intend to prevent 
the accumulation of truncation error. For these purposes 
we employ a functional I — J Cdt, which is not only sta- 



tionary but is also minimum for the actual time evolution 
from ^(tf)) to \ty(t F )). Minimization for the integral 
of the Lagrangian like function £, which is bilinear in 
\^f>(t)), for the time span t r < t < t F draws a way of 
improving MPS on each time slice iteratively. In a sense 
the optimization process that we explain in the following 
can be regarded as 'the finite-time DMRG algorithm', 
which sweeps MPS between t I to t F iteratively. In con- 
trast to the spacial sweeping process in the finite-size 
DMRG algorithm applied to ground state problems, the 
sweeping toward and backward the time direction can be 
performed simultaneously by parallel computation. 

In the next section, we introduce a kind of action / = 
J Cdt that is simply written as the square error with 
respect to the small time evolution by transfer matrices. 
In §3 we explain how the least action principle draws 
the optimization conditions for MPS on each time slices. 
Conclusions are summarized in the last section. 

§2. Least Action Principle 

Consider the time evolution of an isolated quantum 
state in the Shodinger picture 



d 



m)) 



H(t) 



(2.1) 



dt ih 
where H (t) is the Hamiltonian of the system. We have 
divided the standard formulation by ih in order to con- 
centrate on the time evolution of the quantum state. The 
formal solution of this equation from the initial state 
I *(*/)) is g iven by 

!*(**■)> = exp f i fj H(t)dtj Mh)) , (2.2) 

where we assume the time order in the exponential of 
the integral. Normally the function 

d H(ty 



dt 



!*(*)> (2-3) 



is chosen as the Lagrangian. If we do not care about 
the physical dimension of the Lagrangian, there are ac- 
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tually a lot of functions that draws Eq.(2.1) by way of 
the stationary condition. Among such functions, there 
is non-negative one 



£'(*) = (*(*)! (| 



(-- 




(d 


H(t)\ 








a J 



!*(*)> 

(2.4) 



which we treat in the following. Though there is no profit 
of considering £'(t) instead of C in analytical calcula- 
tions, it is of use for finding a variational principle in 
the real-time DMRG method. We define a functional J, 
which corresponds to a kind of action from the initial 
time tj to the final one t F , by the integral 15 ) 



-I'M- 



H(t) 
ih 



!*(*)> 



dt. 



(2.5) 



In order to simplify the formulation, we concentrate on 
time-independent Hamiltonian in the following. (Intro- 
duction of time dependence is straight forward.) 

For the numerical treatment of time evolutions, let 
us divide the time span into N segments and introduce 
discrete time t e = t I + £At, where At = (t F — tj)/N. 
The final state is then formally expressed as 



=exp 



ih 



nh)) = T N mi)) (2-6) 



using the short-time evolution operator 

T = exp(^ff). (2.7) 
In the same manner, the state at t i+1 is written as 

\*(t i+1 ))=TMt i ))=T i + 1 \*(t I )). (2.8) 

As discrete analogues of the functional in Eq.(2.5), we 
employ an error function 



JV-l 



'/ = E 



i=0 



*(t i+1 ))-T|*(t i )) 



or the similar one 

JV-l 



E 

i=0 



T-'mti+iV-mu)) 



(2.9) 



(2.10) 



where we have introduced a backward small-time evolu- 
tion operator 



(2.11) 



These two error functions are actually the same, since 
we have TT _1 = id, and thus 

|*(t i+1 )> -r\^( ti )) = T(r-^{t l+1 )) |*( ti )» 

(2.12) 

is satisfied. Let us focus on the minimization of a part 
of the error 



T-'Mti+iV-Mti)) 



+ 



(2.13) 

that is related to ^(t^)). The stationary condition with 



respect to the variation VE^} 
optimal condition 



l*(*i)>=2 



T- 1 |*(i i+1 )) + T|*(V 1 )) 

for those states at t i < t F , and 

|*(i JV )) = T|*(t JV _ 1 )> 



+ | (5*) draws the 
(2.14) 



(2.15) 



for |*(t F )>. 

As an example of one-dimensional (ID) systems, let us 
consider a lattice model of length L, whose Hamiltonian 
is written as the sum of nearest neighbor interactions 

L-l 

h = E = E + E 

i— £ £— odd £— even 

= H 1 +H 2 , (2.16) 

where we assume the open boundary condition through- 
out this article. In such a case the time evolution is well 
approximated by the Trotter formula 

1 N 



Mt F )) 



exp^^exp^; 



!*(*/)> 



JV 



!*(*/)>. 



(2.17) 



where 7^ and T 2 are written as product of non- 
overlapping local time evolutions 

%. = n ex p(jl- = n T «+i' 

fcodd " £=odd 



% = n ° x p(-^ vn) = n r ^^+i • ( 2 - is ) 

^-even £—cvcn 

The Trotter decomposition introduces a new state 
l^+i/2) = ^1 !*(*<)> between |*.) = |*(t.)) and 
l*i+i> = l*(*i+i)>- The index i + 1/2 of |0. +1 , 2 ) just 
means that it is between i and i + 1; note that |0 i+1 / 2 ) 
does not correspond to |^(^ + At/2)). Now we have 2iV 
numbers of states, and the error function in Eqs.(2.9) or 
(2.10) is modified as 

JV-l 



/ = E |i<w>- T ii*<> 

i=0 
JV-l 

+ E |i^+i)- r 2 1^+1/2) 



i=0 



(2.19) 



As we have done for Eqs.(2.13) and (2.14), minimization 
of this error function draws the following conditions 

1 



T 2 |^_ 1/2 ) + (T, ) 1 1^ +1/2 > 

^ii* i ) + (^r 1 i*i + i> 



2 

'. 1/2) = 2 

for i < N. Only at the final time 

I*JV) = % l0JV-l/2> 



(2.20) 
(2.21) 



should be satisfied. From equation (2.20), it is appar- 
ent that we can deal \4> i+1 / 2 ) equivalently with l^} as 
long as the minimization of the error function / is con- 
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cerned. We therefore explain the optimization of l^) 
only in the following. This variational formulation is 
time-symmetric, in the sense that both 7[ and T 2 are 
invertible and thus I b = If. 

§3. Optimization of Matrix Product States 

In the framework of DMRG method, all the states are 
(implicitly) written as the MPS. Let us write |^) = 
in the form of the orthogonal MPS 3 ~ 5) 



division can be shifted arbitrarily by way of the rela- 

tion 4,5,12,13) 



B[s l }...B{s l ] A B B[ Se+1 }...B[s L ]\ Sl 



(3.1) 



which corresponds to the division of the lattice into the £- 
site left part and the L— I site right part. The notation Sj 
represent the site variable — say, the spin variable — at 
position j. Suppose that the degree of freedom of each Sj 
is d. Since we are treating a system with open boundary 
condition, B\s{[ is a 1 x d matrix (= row vector) 

B xa [ Sl } = 8(a, Sl ) (3.2) 

and B[s L ] is a d x 1 matrix (= column vector) 

B.xkl =S(a,s L ), (3.3) 

where the matrix indices ' x ' and a represent the 1-state 
and rf-state auxiliary variables, respectively. In the fol- 
lowing we consider the case where the size of the matrix is 
less than a fixed integer m, which is the maximum num- 
ber of states kept in the DMRG method. Figure 1 shows 
the structure of the orthogonal MPS. The white circles 
represent spin variables s i , and the black squares repre- 
sent auxiliary variables whose sum is taken over. The 
triangles correspond to the orthogonal matrices B[sj]. 
The large circle at the dividing point of the system rep- 
resent the diagonal matrix A B . 



Fig. 1. Graphical representation of the MPS in Eq. (3.1) for the 
case L = 12. 



The matrices B[sj] in the left side of A B satisfies the 
orthogonal relation 

£ (B af3 [sj])* B ai [ Sj ] = 5 M (j < I) , (3.4) 



and that in the right side of A B satisfies 

^(S a/3 [ Sj ])*B 7/3 [ Sj ]=5 a7 (j>£). (3.5) 

The diagonal matrix A B contains the singular values A B 
with respect to the left-right division of the state I*,) 
at the position I. The matrix A B is dependent to £, and 
satisfies the condition 

Tr(A B ) 2 = ]T(A B ) 2 = l (3.6) 

a 

when |^-) is normalized. The position of the left-right 



B[ S ,]A B = *J S ,]=A B B[ S£ 



(3.7) 



where B[s e ] of B[s e ] A B satisfies Eq.(3.4), and that of 
A B B[sf\ satisfies Eq.(3.5). (See Fig. 2.) Using the renor- 
malized wave function [s e ] defined in the above equa- 
tion, we can express 1^) as 

B[ Sl ] . . . Bls^} * 4 [s e ] B[s e+1 ] . . . B[s L ] \ Sl ...s L ). 

(3.8) 



Fig. 2. Rcnormalized wave function defined in Eq.(3.7). 



Because of the orthogonality in Eqs.(3.4) and (3.5), 
the inner product (^l^) can be expressed simply as 

<*il*i> = £ (^W)'^W 

af)s. 



£ (B al3 [s e ] A B )* (B a/3 [s e ] A B ) 

a/3s e 

£(a^) 2 = i. 





(3.9) 



The right hand side of the first line can be regarded as 
the norm of the renormalized wave function "J^ [s e ] , which 
we express as (^J^J in the following. 

In the same manner as l^), the states 1 0^ 1/2) ana - 

l^i+1/2) can ex P resse d in the form of the orthogonal 
MPS 

l^-i/2) = • • • A[s f ] A A A[s t+1 ] . . . A[s L ] \ Sl ... s L ) 

I W> = C \*i\ ■ ■ ■ C ^ A ° C ^+J ■ ■ ■ ■ ■ • s l) > 

(3.10) 

respectively, and the corresponding rcnormalized wave 



-1/2L 



A[s e ]A A and 



i+i/2l s e\ — 



functions are 
C[s t \ A c . 

Consider a variation of l^) with respect to the local 
change in MPS caused by ^[s e ] -> ^[s e ] + S^[s e }. What 
should be minimized with respect to this variation is 



1/2/ 



l*i> 



(^)" i i^ + i/2)-i^)r+i T 2^-i 

= (*il*i> - (^ + i/2l^i I*,) - to-i/alte ) _1 |*i> + 1 
+ h.c. . (3.11) 

Using the matrix product structures of each state, the 
inner products {4> i+1 / 2 

can be calculated rapidly of the order of m 3 L in compu- 
tational time. 11 ' 14 ) To simplify the notation let us intro- 
duce new renormalized wave functions [s e ] and [s e ] 
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that satisfies the relations 

(*"!**) = E (^wj'^w 

a/3s < 

= <0 i -i/ 2 K^r 1 i*i). 

(* + l*i)=E(*^W)**a/jW 

= <0 i+ i/ 2 l^il*i)- (3-12) 

It is easily shown that the computational cost of obtain- 
ing both fy~[s e ] and ^ + [s £ ] is also of the order of m 3 L. 




h 

Fig. 3. Rcnormalizcd wave functions *J/ — [s^] (upper) and ^[sj 
(lower). Squares represent the local evolution operator t- - +1 
contained in T 2 . 



computation. After improving l"]^) and ^,±1/2) on all 
the time slices alternatively for numbers of times, we ob- 
tain the matrix product states that minimizes the error 
function in (2.19). 16 ) A more realistic way of calculation 
is to boost the state from \^{t N )) to \^{t N+1 )} using the 
conventional real-time DMRG methods, and to optimize 
the obtained MPSs for M(« N) numbers of time slices 
from the frontier after each time boost. 

The explained procedure does not change the size m 
of each matrix. Occasionally it is better to change the 
value of m site by site according to the truncated weights 
in the renormalization process. Variation with respect to 
the extended renormalized wave function 

hlsi'l+i] = 5 M ^B[s e+1 ] = B[s g ]B[s g+1 ] A B 

(3.17) 

makes it possible to adjust m dynamically during the 
calculation. 

Fig. 4. Extended renormalized wave function in Eq.(3.17). 



The quantity in Eq. (3.11) can be then written in short 
- - + l} + h.c. , (3.13) 

and we finally obtain the optimal condition 

* i [s i ) = l*-[s i ] + l* + [s e ] (3.14) 

that improves the renormalized wave function ^[sj. 
Only at the final time t F — t N the condition is modi- 
fied as 

*N[st\ = \*-[*t\- (3-15) 

Even when the states and \4> i±1 / 2 ) are normal- 
ized, the improved ^[s^] created by Eqs.(3.14) and 
(3.15) does not always satisfy the normalization condi- 
tion (^jl^i) = 1. Thus we normalize ^[s^] as 

*i[^]/(*il*i) 1/2 -*i[*/] (3-16) 

after each local optimization. Schmidt orthogonalization 
of thus improved vf^ [s e ] gives improved orthogonal ma- 
trix B[s(\ and the singular value A B at the position I on 
the time slice t i . Sweeping this improving process from 
i = 1 to i = L for several times, one obtains optimized 
with respect to the fixed \<P i±1/2 )- n ' 14) 
To perform such a sweeping process for all the time 
slices requires a huge amount of numerical calculation. 
But this is not totally unrealistic, since improvements 
of I*.) = and 1*^.) = ^ i.)) can be per- 

formed independently under the condition that (^,±1/2) 
for every i is fixed. Also we can say that ^±1/2} can 
be improved simultaneously for every i when all the 
are fixed. The nature enables us to perform the parallel 



§4. Conclusion and discussion 

We reformulate the real-time DMRG method so that 
it minimizes a kind of discrete action, which corresponds 
to the square of error in the numerical time evolution. As 
a result the time symmetry is recovered from the level of 
variational formulation. The minimization process can 
be performed via parallel computation. 

Let us discuss the origin of the slight time asymmetry 
in the conventional real-time DMRG algorithms. Con- 
sider the multiplication of T x to |^) that is expressed 
as the MPS. It is possible to represent T x l*^) exactly in 
the form of MPS again, but this requires more degrees 
of freedom than the original MPS representation of l^). 
Therefore a truncation (= renormalization) process R is 
applied to T x 1^) to keep the degree of freedom nearly 
constant. As a result, very small but non-zero discrep- 
ancy 

TM-RTM (4.1) 

arises. Since the truncation process satisfies R 2 = R, 
and since is obtained by applying R to T 2 , 2 ) , 
we obtain R\^i) = l^). Thus the above discrepancy 
can by written as 

T 1 R\af i )-RT 1 \* i ) = [T 1 ,R]\* i ). (4.2) 

This tiny error cannot be recovered by the multiplication 
of (T 1 ) from the left, since {T x ) RT X is not the iden- 
tity operator. For T 2 there is the same kind of discrep- 
ancy [T 2 , R] \4>i_i/ 2 )- Although these are tiny errors, 
repeated use of the truncation R after each time evo- 
lution may introduce exponentially growing numerical 
error. The variational treatment presented in this arti- 
cle recovers the time symmetry by considering Eq.(3.11). 
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The time boost at t F is still asymmetric as shown in 
Eq.(3.15), and this asymmetry should be removed after- 
ward through the repeated optimization processes. 

In our formulation, we assumed that the Hamiltonian 
can be decomposed into two non-overlapping parts. For 
the cases where there are long-range interactions, the 
multi-targeting scheme by Feiguin and White is of use. 9 -* 
The outline is to optimize ^(i-)) so that the error func- 
tion in Eq.(2.13) is minimized. This process is realized 
by rewriting the involved states ^(^-i)}, "J 7 ^)), and 
\^(t i+1 )) using the same orthogonal matrices. 17 ) It is in- 
teresting that the method of Feiguin and White performs 
the same kind of optimization for the succeeding 4 states 

l*(*i)>> l*(*i+l/3)>. l*(*i+2/3)>' l*^+l))- 18) . 

Remember that there are many possible choices of 
Lagrangian-like function even in continuous time formu- 
lations, and there are much more for the discrete time 
cases. The authors have just considered one of them. 
Looking at the classical simulation of Newton equation, 
there are various techniques that would be implemented 
in real-time DMRG algorithm. For example, automatic 
adjustment of At, use of the symplectic structure, posi- 
tion dependent choice of At, etc. It is worth considering 
what kind of Lagrangian or error function exist toward 
these extensions. 

The variational method introduced in this article 
can be applied for transfer matrix problems in two- 
dimensional classical lattice models, if T -1 is represented 
as a product of local factors. In the same manner, classi- 
cal simulations of quantum circuits 19, 20 - ) are in our scope, 
since the transfer matrices that represent quantum oper- 
ations are always invertible. For those cases where there 
is no inverse, or non-local terms appears in T _1 , we still 
do not obtain an appropriate variational functional. 

T. N. is partially supported by a Grant-in- Aid for Sci- 
entific Research from the Ministry of Education, Science, 
Sports and Culture. 
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15) Note that the integral does not have the proper dimension of 
the action. It is just the functional integral of some function 
C'(t). 

16) After the complete improvement, vp^sj is very close to 
They become the same both in the limit At — > or for suffi- 
ciently large m. 

17) E. Jeckelmann: Phys. Rev. B 66 (2002) 045114. 

18) After expressing the succeeding 3 states by use of the same 
orthogonal matrices, there is no time asymmetry in the time 
boost in Eq. (2.8). However, very few but nonzero error enters 
during the rewriting process of MPS, which plays the same 
role as R in Eqs.(4.1) and (4.2). 
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